% Estimate Bid-Ask illiquidity spread

% A -- is my bid_ask_matrix of bid-ask prices 
% odd columns == bid prices
% even columns == ask prices

% Get the number of columns in A
[rows, cols] = size(A);

% Extract all the odd-numbered columns == Bid prices
bid = A(:, 1:2:cols);

% Extract all the even-numbered columns == Ask prices
ask = A(:, 2:2:cols);

% Estimate bid ask-spread

bid_ask_spread= (ask - bid) ./ ask;  % Calculate the spread using (Ask - Bid) / Ask

% Monthly averages of bid-ask

y = 1;
for i = 0:21.725:2607
BA(y,: ) = mean(bid_ask_spread(i+1:i+21.725,:));  
y=y+1;
end

% Reshape the bid_ask_spread into single column vector 

liquidity=BA(:); 
%_________________________________________________________________________


% Get the log returns of CDS 

r= diff(log(cds));

% Get the Squared returns

R=r.^2; % Squared log returns

% Get the RV of squared log returns for each month (120) 
%_________________________________________________________________________
y = 1;
for i = 0:21.71666:2606
RV_m(y,: ) = sum(R(i+1:i+21.7166,:));  
y=y+1;
end

% Reshape RV into single column vector
RV=RV_m(:);

% Get the Bipower Variation BV sum (from 2 to M)
%__________________________________________________________________________

rr=abs(R);
p1=rr; % rr is the absolute log-return
y=1;
for i=0:1:2606
    C1(y,:)=prod(p1(i+1:i+2,:)); % I need C1 only for to calculate C - drop C1 
    y=y+1;
end

p2=C1;
y=1;

for i = 0:21.71666:2606
C(y,:) = sum(p2(i+2:i+21.71666,:)); % this is the sum (from 2 to M)that I need for BV
y=y+1;
end

clear ('C1', 'p1', 'p2')

mu1= (pi/2); % eq 4 or see Appendix A

BV1 = mu1.^(-2) .* C;
BV=BV1(:);
% TQ = Compute Tripower Quarticity TQ sum (from 3 to M)
%__________________________________________________________________________

p3=rr;
y=1;
for i=0:1:2606
    E1(y,:)=prod(p3(i+1:i+3,:)); % I need E1 only for to calculate E- drop E1 
    y=y+1;
end

p4=E1;
y=1;

for i = 0:21.7:2604
EE(y,:) = sum(p4(i+3:i+21.71666,:)); % this is the sum (from 3 to M)that I need for TQ (TP)
y=y+1;
end

clear ("y", "E1", "p3", "p4", "i")
M = 21.71666;
delta = 1./M;

mu = 1./(4.*delta.*((gamma(7/6))./(gamma(1/2))).^3); % see appendix A

TQ = mu.*EE; % E or EE is the: R = r.^(4/3); the abs(return).^4/3 %see also TP appendix A
TP = TQ(:);
% RJ = Compute Realized Jumps 
%_________________________________________________________________________
 RJ = (RV-BV)./RV;

% Z1 = Compute the Z-test statistics for assessing the significance of the daily
% jump component. Thus, only statistically extreme positive values of RV-BV/RV
% are attributed to the jump component. BNS(2004), Huang and Tauchen(2005)
% and Andersen, Bollerslev and Diebold(2007) for significant jumps test
% statistics based in ratio statistics
%_________________________________________________________________________________________

Z1 = RJ./(((pi/2).^2 + pi - 5)* delta * max(1, TP./(BV.^2) )).^(1/2);

% Realized Jump Volatility == Computes the stat significance of Jumps
%__________________________________________________________________________
% 95% significance level
for i = 1:1:size(Z1,1)
    for j = 1:1:size(Z1,2)
        if Z1(i,j) >  1.644854 %1.644854  is of 95% 3.090232 %99.9%[the correct] % 2.326348 is 99% significance level
RVJ(i,j) = sqrt(RV(i,j))-sqrt(BV(i,j));
RVC(i,j) = sqrt(BV(i,j));
        end
    if Z1(i,j) < 1.644854 
   RVJ(i,j) = 0;
    RVC(i,j) = sqrt(RV(i,j));     
    end
    end
end

% 99 % Significance level
for i = 1:1:size(Z1,1)
    for j = 1:1:size(Z1,2)
        if Z1(i,j) >  2.326348 % 2.326348 is 99% significance level
RVJ_99(i,j) = sqrt(RV(i,j))-sqrt(BV(i,j));
RVC_99(i,j) = sqrt(BV(i,j));
        end
    if Z1(i,j) < 2.326348 
   RVJ_99(i,j) = 0;
    RVC_99(i,j) = sqrt(RV(i,j));     
    end
    end
end

% Compute Co-Jumps Bollorslev, Law and Tauchen (J. of Econometrics 2008)
% % Compute Mean Cross Product (mcp) Statistics (defined as normalized sum of
% individual high-frequency returns for each within day period
%__________________________________________________________________________

n=110; % n is the number of firms on my portfolio
      % M is the number of intraday returns e.g. M=30 
     
mcp_tj = sum(r,2).*(2./(n.*(n-1))); % Eq(11) mcp for each within-day return in the porfolio
mcp=mcp_tj(:);
a=mcp_tj;

y = 1;
for i = 0:21.71666:2606
month(y,: ) = mean(r(i+1:i+21.71666,:)); % monthly average of the mcp 
y=y+1;
end


y = 1;
for i = 0:21.71666:2606
smcp_tj(y,: ) = std(a(i+1:i+21.71666,:)); % 21.71666:2606 monthly standard deviation of the daily mcp 
y=y+1;
end

b=repmat(month',1,22);
c=b(:); %delete last n observations to match the length with mcp

d=repmat(smcp_tj',1,22);
e=d(:); %delete last n observations to match the length with mcp

av_mcp_tj_1=repmat(month',151,1);
av_mcp = reshape(av_mcp_tj_1,[],1);

smcp_tj_1=repmat(smcp_tj',151,1);
smcp = reshape(smcp_tj_1,[],1);

% Compute the z-stats for co-jumps Eq (14) daily/weekly/monthly
% _________________________________________________________________________
zmcp_tj = (mcp_tj - c)./e;

y = 1;
for i = 0:21.71666:2606
month_zmcp(y,: ) = mean(zmcp_tj(i+1:i+21.71666,:)); % monthly average of the mcp 
y=y+1;
end

zfirms=repmat(month_zmcp',1,110);

z=zfirms(:);
